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Abstract. We present and analyze stochastic nonlinear differential equations 
generating signals with the power-law distributions of the signal intensity, 1//^ noise, 
power-law autocorrelations and second order structural (height-height correlation) 
functions. Analytical expressions for such characteristics are derived and the 
comparison with numerical calculations is presented. The numerical calculations 
reveal links between the proposed model and models where signals consist of bursts 
characterized by the power-law distributions of burst size, burst duration and the inter- 
burst time, as in a case of avalanches in self-organized critical (SOC) models and the 
extreme event return times in long-term memory processes. The presented approach 
may be useful for modeling the long-range scaled processes exhibiting 1// noise and 
power-law distributions. 
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1. Introduction 

The inverse power-law distributions, autocorrelations and spectra of the signals, 
including 1/f noise (also known as 1/f fluctuations, flicker noise and pink noise), as 
well as scaling behavior in general, are ubiquitous in physics and in many other fields, 
counting natural phenomena, spatial repartition of faults in geology, human activities 
such as traffic in computer networks and financial markets. This subject is a hot research 

topic for many decades (see, e.g., [HEIEIIUEIIHIITIIEIIHI ^ anc ^ re f erences herein). 
An up-to-date bibliographic list on 1/f noise of more than 1300 papers is composed by 
Wentian Li [TT| . 

Widespread occurrence of signals exhibiting such a behavior suggests that a generic, 
at least mathematical explanation of the power-law distributions might exist. Note that 
the origins of two popular noises, i.e., the white noise - no correlation in time, the power 
spectrum S(f) ~ f°, and integral of the white noise - the Brownian noise (Wiener 
process), no correlation between increments, the power spectrum S(f) ~ f~ 2 , are very 
well known and understood. 1/f 13 noise with < < 2, however, cannot be realized 
and explained in a similar manner and, therefore, no generally recognized explanation 
of the ubiquity of 1/f noise is still proposed. 

Despite numerous models and theories proposed since its discovery more than 80 
years ago p], the intrinsic origin of 1/f noise still remains an open question. Although 
in recent years it is annually published about 100 papers with the phrases '1// noise', 
'1// fluctuations' or 'flicker noise' in the title, there is no conventional picture of the 
phenomenon and the mechanisms leading to 1/f fluctuations are not often clear. Most 
of the models and theories have restricted validity because of the assumptions specific 
to the problem under consideration. Categorization and summary of the contemporary 
stage of theories and models of 1/f noise are rather problematic: on one hand, due 
to the abundance and variety of the proposed approaches, and on the other hand, for 
the absence of the recent comprehensive review of the wide-ranging "problem of 1/f 
noise" and because of the lack of a survey summarizing the current theories and models 
of 1/f noise. We can cite only a pedagogical review of 1/f noise subject by Milotti 
[12] . Wentian Li presents some kind of classification by the categories of publications 
related with 1/f noise until 2007 [13]. In the peer- reviewed encyclopedia Scholarpedia, 
P3] there is also a short current review on the subject under the consideration. Thus, 
we present here only a short and partial categorization of 1/f noise models with the 
restricted list of references. 

Until recently, probably the most general and common models, theories and 
explanations of 1/f noise have been based on some formal mathematical description 
such as fractional Brownian motion, half-integral of the white noise or some algorithms 
for generation of signals with the scaled properties [15] and popular modeling of 1/f 
noise as the superposition of independent elementary processes with the Lorentzian 
spectra and proper distribution of relaxation times, e.g., l/r re i ax distribution [IB] . The 
weakness of the later approach is that the simulation of 1/f 13 noise with the desirable 
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slope (3 requires finding the special distributions of parameters of the system under 
consideration, at least a wide range of relaxation time constant should be assumed in 
order to correlate with the experiments (see also El UHl [T7]). 

Models of 1/f noise in some particular systems are usually specific and do not 
explain the omnipresence of processes with 1//^ spectrum. Predominantly 1/f noise 
problem has been analyzed in conducting media, semiconductors, metals, electronic 
devices and other electronic systems [Tj (5J [161 EH] ■ The topic of 1/f noise in such 
systems has been comprehensively reviewed [5], even recently [6]. Nevertheless, despite 
numerous suggested models, the origin of flicker noise even there still remains an open 
issue: "More and more studies suggest that if there is a common regime for the low 
frequency noise, it must be mathematical rather than the physical one" [6]. Here we 
can additionaly mention the disputed quantum theory [19] of 1/f noise and satisfactorily 
interpreted 1/f noise in quantum chaos [20] . 

In 1987 Bak, Tang, and Wiesenfeld [211 122] introduced the notion of self-organized 
criticality (SOC) with one of the main motivation to explain the universality of 
1/f noise. SOC systems are nonequilibrium systems driven by their own dynamics 
to a self-organization. Fluctuations around this state, the so-called avalanches, are 
characterized by the power-law distributions in time and space, criticality implying 
long-range correlations. The distributions of avalanche sizes, durations and energies are 
all seen to be power laws. 

Two types of correlations should be distinguished in SOC: the scale-free distribution 
of their avalanche sizes and temporal correlations between avalanches, bursts or (rare, 
extreme) events. In the standard SOC models the search of 1/f 13 noise is based on the 
observable power-law dependence of the burst size as a function of the burst duration 
and the power-law distribution of the burst sizes, with the Poisson distributed interevent 
times. Such power-laws usually result in the relatively high-frequency power-law, 1/f^, 
behavior of the power spectrum with the exponent 1.4 < < 2 [23~1 124"] . This mechanism 
of the power- law spectrum is related to the statistical models of 1/f^ noise representing 
signals as consisting of different random pulses [251 EEl 127] . 

It should also be mentioned that originally SOC has been suggested as an 
explanation of the occurrence of 1/f noise and fractal pattern formation in the dynamical 
evolution of certain systems. However, recent research has revealed that the connection 
between these and SOC is rather loose [28]. Though an explanation of 1/f noise was 
one of the main motivations for the initial proposal of SOC, time dependent properties 
of self-organized critical systems had not been studied much theoretically so far [29] . 

It is of interest to note, that paper |21j is the most cited paper in the field of 1/f 
noise problem, but it has been shown later on [231 121] that the proposed in Ref [2TJ 
mechanism in SOC systems results in 1/f 13 fluctuations with 1.5 < f3 < 2 and does not 
explain the omnipresence of 1/f noise. On the other hand, we can point a recent paper 
[30] where an example of 1/f noise in the classical sandpile model has been provided. 

It should be emphasized, however, that another mechanism of 1/f 13 noise, based 
on the temporal correlations between avalanches, bursts or (rare, extreme) events, may 
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be the source of the power-law 1/f 13 spectra with (3 < 1 [31]. Moreover, SOC is closely 
related with the observable 1/f 13 crackling noise [32], Barkhausen noise [33], fluctuations 
of the long-term correlated seismic events |34j and 1/f 13 fluctuations at non-equilibrium 
phase transitions [35] . 

Ten years ago we proposed [EJ E], analyzed [36] and later on generalized pH] a 
simple point process model of 1/f 13 noise and applied it for the financial systems [37J. 
Moreover, starting from the point process model we derived the stochastic nonlinear 
differential equations, i.e., the general Langevin equations with a multiplicative noise 
for the signal intensity exhibiting 1/f 13 noise (with 0.5 < (3 < 2) in any desirable wide 
range of frequency / [38] . Here we analyze the scaling properties of the signals generated 
by the particular stochastic differential equations. We obtain and analyze the power-law 
dependencies of the signal intensity, power spectrum, autocorrelation functions and the 
second order structural functions. The comparison with the numerical simulations is 
presented. 

Moreover, the numerical analysis reveals the second (reminder, that we start from 
the point process) structure of the signal composed of peaks, bursts, clusters of the events 
with the power-law distributed burst sizes S, burst durations T and the inter-burst time 
9, while the burst sizes are approximately proportional to the squared durations of the 
bursts, S ~ T 2 . Therefore, the proposed nonlinear stochastic model may simulate SOC 
and other similar systems where the processes consist of avalanches, bursts or clustering 
of the extreme events [211 [22l [231 EH EH [291 EHl EU [321 [331 EH [35l [39]. 



2. The model 



We start from the point process, 

x(t) = a^25(t-t k ), (1) 

k 

representing the signal, current, flow, etc, x(t), as a sequence of correlated pulses or 
series of events. Here S(t) is the Dirac 5-function and a is a contribution to the signal 
x(t) of one pulse at the time moment t k . Our model is based on the generic multiplicative 
process for the interevent time r k = t k+ i — t k , 

Tfc+i = r k + 7rf' t ~ 1 + ar^e k , (2) 

generating the power-law distributed, 

P k (T k )~ T «, « = ^-2/i, (3) 

sequence of the interevent times r k [TU1 [37]. 

Some motivations for equation (T5|) were given on papers [8j [TQl |36l [37] . Additional 
comments are presented below, after equation 

Therefore, in our model the (average) interevent time r k fluctuates due to the 
random perturbations by a sequence of uncorrelated normally distributed random 
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variables {e^} with zero expectation and unit variance, a is the standard deviation 
of the white noise and 7 1 is a coefficient of the nonlinear damping. 

Transition from the occurrence number k to the actual time t in equation (j2J) 
according to the relation dt = r^dk yields the Ito stochastic differential equation for the 
variable r(t) as a function of the actual time, 

dr = 7 r 2 ^ 2 dt + ar^- 1/2 dW, (4) 

where W is a standard Wiener process. Equation (j3J) generates the stochastic variable 
r, power-law distributed, 

^)^r~^ (5) 
in) 

in the actual time t. Here (rfc) is the average interevent time. r(t) may be interpreted as 
the average time- dependent interevent time of the modulated Poisson-like process with 
the distribution of the interevent time 

Ppfo) = -L e-T,/T(t) = n (t) e~ n ^, (6) 

T{t) 



where n{t) = l/r(t) is the time dependent rate of the process 

Additional support for the stochastic model <^ - (E]) of the scaled processes and 
1/fP noise is the following. The fluctuations of the intensity of the signals, currents, 
flows, etc, consisting of the discrete objects (electrons, photons, packets, vehicles, 
pulses, events, etc) are primarily and basically defined by fluctuations of the (average) 
interevent, interpulse, interarrival, recurrence, or waiting time. Equation (jjj) is a special 
case of the general non-linear Langevin equation 

dr = d(f)dt + b(r)dW(t) (7) 

with the drift coefficient d(r) and a multiplicative noise &(r)£(t) for the (average) 
interevent time r(t), with £(t) being a white nose defined from the relation dW(t) = 
£(£)d(t). Equation (J7|) is a straight analogy of the well-known Langevin equation for 
the continuous random variable x. For the process consisting of the discrete objects the 
intensity of the signal fluctuates due to fluctuations of the rate, i.e., density of the objects 
in the time axis, which is a consequence of fluctuations of the interarrival or interevent 
time. Equation ([7]) in reality represents (in the simplest form) such fluctuations due to 
random perturbations by white noise. 

In papers [H HI [10] it has been shown that the small interevent times and clustering 
of the events make the greatest contribution to noise, low frequency fluctuations 
and exhibition of the long-range scaled features. Therefore, it is straight to approximate 
the non-linear diffusion coefficient b(r) and the distribution of the interevent time P*(t) 
in some interval of small interevent times r by the power-law dependences or expensions, 

b(T) = aT lt -i, (8) 
Pt(r) ~ r a+1 . (9) 
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The power-law distribution of the interevent, recurrence, or waiting time is 
observable in different systems from physics and seismology to the Internet, financial 
markets and neural spikes (see, e.g., pHl ETJ |39]). It should be noted that the 
multiplicative equations with the drift coefficient d{r) proportional to the Stratonovich 
correction for the drift, leading the transformation from the Stratonovich to the It 6 
stochastic differential equation [ID], i.e., when 

d(r) ~ h(r)b'(r), (10) 

with the power-law depending, like ((HD, diffusion coefficient b(r), generates the power- 
law distribution of the stochastic variable. Equations (T5J) and (jl]) are definitely of such 
kind. Therefore, equation (JI]) is one of the simplest multiplicative equations for the 
interevent time, modeling scaled processes, while equation (2) is just the lowest order 
difference equation following from equation (jlj) when the step of integration Atk equals 
the interevent time 

The Ito transformation in equation (jlj) of the variable from r to the averaged over 
the time interval r intensity of the signal x(t) = a/r(t) [3H] yields the class of Ito 
stochastic differential equations 

dx = (n - l-X) x 2r, ~ l dt s + x"dW (11) 



2 

for the signal as function of the scaled time 



° 2 , ^ 

-r-o-t. 12 



Here the new parameters 



r/=--/i, A = 3 + a = 4 + 2(77 — 1) (13) 



27 

2 ' ' ' " a 2 

have been introduced. 

The Fokker-Plank equation associated with equation (TTTT) gives the power-law 
distribution density of the signal intensity 

p(x) ~ 4c ( m ) 

x x 

with the exponent A. 

For A > 1 distribution ( j!4p diverge as x — > 0, and, therefore, the diffusion of x 
should be restricted at least from the side of small values, or equation (ITT]) should be 
modified. Thus, further we will consider the modified equation for x > only, 

dx=L- ^ (x m + xf^dts + (x m + x) r 'dW, (15) 

with the additional small parameter x m restricting the divergence of the power-law 
distribution of x at x = 0. 

Equation ( fl5l) for small x x m represents the linear additive stochastic process 
generating the Brownian motion with the steady drift, while for x 3> x m it reduces to 
the multiplicative equation (ITT]) . 
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3. Analysis of the model 



P(x) = | T A , x > 0, A > 1. (16) 



The Fokker-Plank equation associated with equation ( 1151) gives the steady-state solution 
for distribution of x, 

( A - l)^ - 1 

{Xm ~\~ X) 

We can obtain the power spectral density of the signal generated by equation (JTBI) from 
equation (28) derived in paper [10J. After some algebra we can write 

s{f) = jj, />>A = ^T^^ 1) (17) 



with 



A (A - l)r(/3 - 1/2)**- 1 /2 + A - 2^ ; - 



and 



20r(^- l)sin(7r/3/2) V 2vr 

for 0.5 </3<2, 4 — 77 < A < 1 H- 2?7 and r\ > 1. Note that the frequency / in equation 
( ITT)) is the scaled frequency matching the scaled time t s (fT2)) . 

The autocorrelation function C(s) of the process can be expressed according to 
Wiener-Khinchin theorem as the inverse Fourier transform of the power spectrum, 

POO 

C(s) = (x(t)x(t + s))= / S(f) cos(27r/ S )d/. (20) 

A pure 1/ power spectrum is physically impossible because of the total power 
would be infinity Depending on whether is greater or less than one it is necessary 
to introduce a low frequency cutoff f m \ n or a high frequency cutoff / max |JTJ H2]. For 
calculation of the autocorrelation function according to equation (12TJ)) when (3 > it is 
not necessary to introduce the hight frequency cutoff. 

Usually one introduces a discontinuous transition to the flat spectrum at the lower 
cutoff /mi n [HI H2] . Here at low frequencies we will insert the smooth transition to the 
flat spectrum in the vicinity of f ~ /mm, i- e -, we w ih approximate the power spectrum 
(H7J as 

8U)=7J *J*!a: (21) 



(/0 2 + r-y /2 \ a/ ?>, / » /„. 

Inserting (|2T)) into equation (|20)) we obtain 

C{*)= *f rB _ 1 &K h {z), (22) 



where Kh{z) is the modified Bessel function, z = 2nfos and h = ((3 — l)/2. 
The first two terms of expansion of equation (1221 in powers of z are 



ar(/?/2)/, 







(23) 
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for h 7^ 0, i.e, 0^1, and 

C(s) = Air (27r/ s) ~ A[-7 - ln(vr/ s)] = const - Alas (24) 

for h = 0, i.e., for the pure 1// noise with (3 = 1. Here 7 ~ 0.577216 is Euler's constant. 

The leading terms of expression (|23|) are different, depending on whether < 1 or 
> 1. Thus for /i < 0, i.e., when < /3 < 1 

c(s) = 727^ sm UJ~^' (25) 

while for h > 0, i.e., for 1 < /3 < 3 

C(s) = C(0) - fls"" 1 . (26) 

Here 

C(0) = <»»)= rS(/)d/= *ffff> (27) 

and 

5 = 2r (|) r ff ;L = - A ^- lT{1 - m si ° (t) ■ (28) 

For = 2 equations (1221 and (j261) - pg|l yield 

/47T 

= C{0)e- 2nfos = ^e~ 27r/ ° s = C(0) - Att 2 s ± ... (29) 
2jo 



It should be noted, that particular cases f[2"4"|) - f[2"9l of the general expressions 
and (123]) are in agreement with the results of papers [4TJ H21 03] obtained with the non- 
uniform cutoff of the spectrum at low frequency. On the other hand, the introduced 
parameter h for > 1 coincides with the Hurst exponent if [42J, 



0, O<0<1, 

H= I 1(0-1), l<0<3, (30) 

1, 3 < < 4. 

The exponent if is associated with the scaling of the second order structural 
function, or height-height correlation function j4TJ 02J H3j E] 

F(s) = F 2 (s) = (\x(t + s)- x(t)\ 2 ) ~ s 2H . (31) 

The exponent H characterizes the power-law diffusion rate, as well. This variance of 
the differenced time series (delayed signal) may be expressed as 

F(s) = (x 2 (t + s)> + (x 2 (t)) - 2 (x{t + s)x{t)) = 2[C(0) - C{s)\. (32) 

Substituting expressions ff20|) and (T27j) into (1321 we have 

POO 

F(s) = F 2 (s) = 4 / S(f) sm 2 (nsf)df. (33) 
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For the convergence of integral in ( 13 3 p at < 1 we need to cut off the power-law 
spectrum fflTl) at high frequency / max - Then the leading terms of the height-height 
correlation function (1331) are 



1-/3 
max 



F(s) 



2Ax( 1-/3 



r(2-/3) 



(2vr/ n 



sm 



7T/3 

~2~ 



< < 1, 



(34) 



[ln(vr/ max s) -7], /3 = 1. 
For 1 < /3 < 3 the integral in (I3"3"|) may be integrated exactly and we have 



-2,41/(1 - )')sm ( ^ ) (2tts)^-' 



2 J Bs /3 " 1 , 1 < p < 3. (35) 



The spatial power spectrum and the height-height correlation function (I3~T]) are 
used for analysis of rough self-affine surfaces and assessing the growth mechanism of 
thin films [45J H6j H7J HH H9], as well. There sometimes the violation of the scaling 
relation = 2H + 1 is observable [4H1 1471 15U] . 



10 3 
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Figure 1. Examples of the numerically computed signals according to equations (|37|) 
and ([38jl with the parameters 77 = 2, x m = and different values of (3 and A: 

/3 = 1/2 when A = 2, /3 = 1 when A = 3 and (3 = 3/2 when A = 4 and the interburst 
time 9j as a function of the occurrence number j of the events peaking above the 
threshold value Xth — 0.1 for the pure 1/f noise with f3 = 1 when A = 3. 



4. Numerical analysis 



For the numerical analysis we have to solve equation (Tl5|) and analyze the obtained 
numerical solutions. We can solve equation ( TT5l) using the method of discretization with 
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10~ 2 10" 1 10° I0 1 10 2 10 3 10~ 2 10" 1 10° 10 1 10 2 10 3 



s s 

Figure 2. Distribution density, P(x), power spectral density, S(f), autocorrelation 
function, C(s), and the second order structural function, F2 (s), for solutions of 
equation (fT5j) with 77 = 2, x m = and different values of the parameter A: A = 2 
(circles), A = 2.5 (squares) and A = 3 (triangles) in comparison with the analytical 
results (solid lines) according to equations ([15]). ([15]) - ([2"2"]l . ([M| and ([55]) . respectively. 




10" 2 10"' 10° I0 1 10 2 10 3 10" 2 10" 1 10" 10 1 10 2 10 3 



s s 

Figure 3. As in figure [2] but for the parameters A = 3 (circles), A = 3.5 (squares) and 
A = 4 (triangles). 
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the variable step of integration 

k 2 

hi = At, = - r (36) 

where k is a small parameter while the exponent I rules the dependence of the integration 
step on the value of the variable x. Thus, I = corresponds to the fixed step, for / = 1 
we have analogy with equation (T5]) when the step is proportional to the interevent time 
Tfc, I = 2 (rj — 1) matches the case when the change of the variable x in one step is 
proportional to the value of the variable at time of the step [38J and so on. As a result 
we have the system of the difference equations 

x i+1 = Xi + k 2 yr] - (x m + Xi) 2 ^ 1 ' 1 + k (x m + Xi) v ~ l/2 Si, (37) 

K 2 

k+i =U + - - v Xi > 0. (38) 

Here Si is a set of uncorrelated normally distributed random variables with zero 
expectation and unit variance. In the Milstein approximation equation (|3"7|) should 
be replaced by the equation 

K 2 K, 2 

x i+ i = Xi + — (rj - A) (x m + Xi) 211 ' 1 ' 1 + k ( m + Xi) v ~ l/2 ^i + ^Y (x m + Xi) 211 ' 1 ' 1 e 2 . (39) 

Numerical analysis indicate that the variable of equation ( jl5j) exhibits some peaks, 
bursts or extreme events, corresponding to the large deviations of the variable from 
the appropriate average value. As examples, in figure [1] we show the illustrations of the 
signals generated according to equation ( TT51) for different slopes of the signal distributions 
and the dependence of the interburst time 6j on the burst occurrence number j. We see 
that the the computed signal is composed of bursts of different size with a wide range 
distribution of the inerburst time. In figures [2] and [3] the numerical calculations of the 
distribution density, P (x), power spectral density, S (/), autocorrelation function, C(s), 
and the second order structural function, F 2 (s) = ^JF (s), for solutions of equation 
(j!5p with rj = 2, x m = 10~ 2 and different values of the parameter A are presented. 
We see rather good agreement between the numerical calculations and the analytical 
results except for the structural function F 2 (s) when A > 3. Numerical evaluation of 
the structural function in a case of the steep power-law distribution is problematic, 
because in the calculation one needs to average (squared) small difference of the rare 
large fluctuations. 

In figure H] we demonstrate numerically that the size of the generated bursts S 
is approximately proportional to the squared burst duration T, i.e., S oc T 2 , and 
asymptotically approximately power-law distributions of the burst size, P(S) ~ S 1-13 , 
burst duration, P(T) ~ J 1 - 1 - 5 an( } interburst time, P{9) ~ #~ L5 , for the peaks above 
the threshold value x t h of the variable x(t). 

It should be noted that the parameter rj = 2 yields in equation flU) the additive 
noise and the linear relaxation of the signal x = a/r, i.e., the simple (pure) Brownian 
motion in the actual time of the interevent time with the linear relaxation of the signal. 
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Figure 4. Dependence of the burst size S as a function of the burst duration T: traces 
from the top of the figure for size 5 above the threshold value x t h — 0.02, 0.1 and 0.5, 
respectively, distributions of the burst size, P(S), burst duration, P(T), and interburst 
time, P(9), for the peaks above the the threshold value Xth =0.1. Calculations are as 
in figures [5] and [3] with the parameters r\ — 2, x m — 10 -2 and different values of the 
parameter A: A = 2 (circles), A = 3 (squares) and A = 4 (triangles). 



5. Conclusions 

Starting from the multiplicative point process we obtain the stochastic nonlinear 
differential equations, which generate signals with the power-law statistics, including 
1//^ fluctuations. We derive analytical expressions for the probability density of the 
signal, for the power spectral density, the autocorrelation function, the second oder 
structural function and demonstrate that the analytical results are in agreement with 
the results of numerical simulations. The numerical analysis of the equations reveals 
the secondary structure of the signal composed of peaks or bursts, corresponding to the 
large deviations of the variable x from the proper average fluctuations. The burst sizes 
are approximately proportional to the squared duration of the burst. According to the 
theory pU E7] such dependence for the uncorrelated bursts should result in 1/ f 13 noise 
with (3 ~ 2 in the relatively high-frequency region. The power-law distribution P{9) 
of the interburst time 9 indicates in correlation of the burst occurrence times and may 
result in 1/ f 13 noise with (3 < 2, similarly to the point process model [HI El HO]. On 
the other hand, the proposed model reproduces 1/ / noise models and the processes not 
only in SOC and crackling systems but it is related with the clustering Poisson process 
[5T~j . 1/f noise due to diffusion of defects or impurity centers in semiconductors |52j, 
1/f noise in nanochannels, single-channel and ion channel currents [53], etc. Therefore, 
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the presented and analyzed model may be used for simulation the long-range scaled 
processes exhibiting 1/f noise, power-law distributions and self-organization. 
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